Unidentified case

INIT

#required packages to make a package
pacman::p_load(devtools, roxygen2, testthat)
#dependecies
#pacman::p_load(tidyverse, data.table, DT, plotly, gridExtra, latex2exp, here, logr, future.apply, furrr)
#pacman::p_load(dplyr,parallel,sandwich,here,future,logr,profvis,lobstr,parallel)
#https://ourcodingclub.github.io/tutorials/writing-r-package/
#create_package
#https://combine-australia.github.io/r-pkg-dev/functions.html
#load_all()
#run these two lines and then library(GridMaker) in normal rstudio to install package
# load_all loads a package. It roughly simulates what happens when a package is installed and loaded with library().
devtools::load_all()

devtools::install()

library(vars)

Loading required package: MASS

Loading required package: strucchange

Loading required package: zoo

Attaching package: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

Loading required package: sandwich

Loading required package: urca

Loading required package: lmtest

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ stringr::boundary() masks strucchange::boundary()
✖ dplyr::filter()     masks stats::filter()
✖ dplyr::lag()        masks stats::lag()
✖ dplyr::select()     masks MASS::select()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(plotly)

Attaching package: 'plotly'

The following object is masked from 'package:ggplot2':

    last_plot

The following object is masked from 'package:MASS':

    select

The following object is masked from 'package:stats':

    filter

The following object is masked from 'package:graphics':

    layout

library(iStabi)

Warning: replacing previous import 'purrr::transpose' by
'data.table::transpose' when loading 'iStabi'

Warning: replacing previous import 'data.table::first' by 'dplyr::first' when
loading 'iStabi'

Warning: replacing previous import 'data.table::last' by 'dplyr::last' when
loading 'iStabi'

Warning: replacing previous import 'data.table::between' by 'dplyr::between'
when loading 'iStabi'

sam <- iStabi::data_sim_biVAR

THE MODEL

True parameters are

\[ \begin{align} y_{1,t} &= \textcolor{blue}{\beta} \,\, y_{2,t} + 0.5 y_{1, t-1} + 0.2 y_{2, t-1} + w_{1, t} \\ y_{2,t} &= \textcolor{blue}{\alpha} \,\, y_{1,t} + 0.1 y_{1, t-1} + 0.4 y_{2, t-1} + w_{2, t} \end{align} \]

\[ \begin{bmatrix}1 & -0.15 \\ 0.75 & 1\end{bmatrix} y_t = \begin{bmatrix}0.5 & 0.2 \\ 0.1 & 0.4\end{bmatrix} y_{t-1} + w_t \]

with \(\beta = 0.15\) and \(\alpha = -0.75\).

ESTIMATION VAR

varModel <- VAR(cbind('i' = sam$y1, 's' = sam$y2), p = 1, type = "const")
Res = residuals(varModel)

# Use for computing the moment function
Omega <- cbind('ii' = Res[, 'i']*Res[, 'i'], 'is' = Res[, 'i']*Res[, 's'],
               'ss' = Res[, 's']*Res[, 's'])

Y <- cbind( 'omega_ii' = Omega[, 'ii'], 'omega_ss' = Omega[,'ss'], 'omega_is' = Omega[,'is'] )

# Compute the moment function at each t. 
# INPUT: a named vector with values for the structural parameters
# OUTPUT: a TxG matrix where each row corresponds to the value of the moment function on that obseration (date)
# an each column the value fot ehfat moment equation
GET_MOM_VALUES <- function(t0, ...) {
  list2env(as.list(t0),envir = environment())
  
  #Y <- cbind( 'omega_ii' = Res[,'i']^2, 'omega_ss' = Res[,'s']^2, 'omega_is' = Res[,'i'] * Res[,'s'] )
  
  #Y %*% c(-theta0['alpha'],
  #        -theta0['beta'],
  #        1 + theta0['alpha'] * theta0['beta'])
  
  #fT 
  as.matrix(-alpha*Omega[, 'ii'] - beta*Omega[, 'ss'] + (1+alpha*beta) * (Omega[, 'is']))
  # fT - lambda
  
  #matrix(fT - mean(fT),
  #       ncol = 1) # (T-varLeng) X 1
  
} 

GET_B <- function(theta0) {
  c(-theta0['alpha'],
    -theta0['beta'],
    1 + theta0['alpha'] * theta0['beta'])
}

SEARCH

Search settings

#### Search settings ----
NCORES = params$NCORES
STEPSPERCORE = 200
# Size of adjacent cells to explores is NCORES * STEPSPERCORE

# PARAMS_CFG a tibble where each row defines properties of the parameters in the name column
# stepSize: size of the step for each parameter. You can have different sizes (and precision) for each paramter.
# minVal and maxVal: Limits of the parameter space for each parameters. A list where each element defines the minimum and maximum for each parameter
# ORDER: the order of the parameter has to match the order in theta (as used by the get_b() fcn)
PARAMS_CFG <- tribble(~name,   ~iniVal, ~stepSize, ~minVal, ~maxVal,
                      'alpha',   0.3,    0.05,        -1,      1,
                       'beta',  -0.2,    0.05,        -1,       1) %>%
  mutate('N_TICKS' = ((maxVal-minVal)/stepSize)+1)

NEQS = 1
MODELEQS <- list(list('Y' = Y,  'get_b'=NULL))
                 #'X'=X, 'Z'=Z))
NCORES = 4
# Level of condfidence of the S set. The closer to 1, the more robust the search is to local minima but it takes longer.
SPVALUE = 0.95
# R is NEQS*ncol(Z) where ncol(Z) is the number of instruments
# By default this matrix has to be automatically built depending on the dimensions of MODELEQS.
# The user can alternatively provide this matrix
R = 1
if (exists('Z', mode='numeric')) {
  R <-  rbind(
        cbind(matrix(1, nrow=ncol(Z), ncol=ncol(Z)))
  )
}

theta0 = PARAMS_CFG$iniVal
names(theta0) = PARAMS_CFG$name
head(GET_MOM_VALUES(theta0))

         [,1]
1 -0.56332302
2  0.16689846
3 -0.04654581
4 -0.05713618
5 -0.18004123
6 -0.02825342

verify_global(PARAMS_CFG, Y)
Times = NROW(Res) # time periods

generate_global(PARAMS_CFG,  SPVALUE)

Will use 4 cores for parallel processing

Global search grid

gridLines <- build_lines(PARAMS_CFG, 'iniVal')

# paramsCfg <- PARAMS_CFG
# thisGridLines <- thisGridLines
# NEQS = 1
# psUnderNull = 1
Grid <- build_grid(PARAMS_CFG, gridLines)

[1] "The grid should have: 1600 rows"

Log is not open.

Test

theta0 = c('alpha' = -0.75, 'beta' = 0.15) # the true theta with the estimated (nuisance) parameter 0.378
comp_tests(theta0, typeVarF='NeweyWest', VffLags = 1)

        S   qLLStab    c2_hat 
0.8425139 3.6393726        NA 

comp_tests(theta0, typeVarF='NeweyWest', VffLags = NULL)

        S   qLLStab    c2_hat 
0.9367893 4.0466100        NA 

comp_tests(theta0, typeVarF='const', VffLags = NULL)

        S   qLLStab    c2_hat 
0.7055631 3.0477918        NA 

comp_tests(theta0, typeVarF='HAC', VffLags = NULL)

        S   qLLStab    c2_hat 
0.9892098 4.2730485        NA 

Start search

#logFileName <- 'log_file.txt'
#log_open(logFileName, logdir=FALSE, compact=TRUE)
#log_print(paste0("Local search? ", params$LOCALSEARCH))

#log_print(paste0("IniRows size (in bytes): ", lobstr::obj_size(IniRows)))
# rover_Perseverance will update the Grid in the Global env
#plan(multisession, workers=params$NCORES)

#log_print("Calling garbage collector before Stage2_rover_Perseverance")
#expResult <- Stage2_rover_Perseverance(iniVal=PARAMS_CFG$iniVal, WTLags = 1, localSearch = FALSE) # updates the Grid by in place
Grid <- readRDS( 'Grid_Unidentified_model.rds')
#  saveRDS(Grid, 'Grid_Unidentified_model.rds')
NFIXPARAMS=1

head(Grid)

Key: <alpha, beta>
   alpha  beta        S  qLLStab ROW_NUMBER c2_hat genS_qLL
   <num> <num>    <num>    <num>      <int>  <num>    <num>
1:    -1 -1.00 254.7320 11.16862          1     NA 242.7432
2:    -1 -0.95 249.3581 11.29589          2     NA 237.9850
3:    -1 -0.90 242.5221 11.36426          3     NA 231.8389
4:    -1 -0.85 234.3027 11.37068          4     NA 224.3731
5:    -1 -0.80 224.8306 11.31425          5     NA 215.7057
6:    -1 -0.75 214.2824 11.19630          6     NA 205.9985

Grid[order(qLLStab)][1:10]

    alpha  beta         S   qLLStab ROW_NUMBER c2_hat  genS_qLL
    <num> <num>     <num>     <num>      <int>  <num>     <num>
 1:  0.05 -1.00 22.836937 0.9207345        841     NA 21.681586
 2:  0.00 -1.00 29.695267 0.9310352        801     NA 27.926732
 3:  0.10 -1.00 17.258594 0.9463928        881     NA 16.636024
 4: -0.05 -1.00 38.028594 0.9873345        761     NA 35.558784
 5:  0.15 -1.00 12.775013 0.9995450        921     NA 12.613193
 6: -0.05 -0.95 22.475247 1.0509994        762     NA 21.483042
 7: -0.10 -0.95 30.765062 1.0541015        722     NA 29.022340
 8:  0.20 -1.00  9.217574 1.0732271        961     NA  9.452839
 9:  0.00 -0.95 15.953995 1.0996115        802     NA 15.603243
10: -0.10 -1.00 48.029135 1.1011542        721     NA 44.764004

BUILD SETS

Critical values

Critical values for \(qLL-\tilde{S}\) for 1 moment condition from table VII Suplementary material, p3 Magnusson Mavroeidis (2014).

THe S-test follows a Chi-Square distribution with 1 degrees of freedom (equatl to the number of moment conditions).

Grid has 1600 rows.

The number of fixed identified parameters is 1.

The number of estimated parameters under the null is 0.

NSIPARAMS = 0 # not needed in the new version of the package
PCONFS = c(0.85, 0.90)
get_CVs <- function() {
  
  CVs <- list('qLLStab'=NA, 'S'=NA, 'genS_qLL'=NA)

  CVs$qLLStab <- filter(iStabi::data_CVs$qLLStab, DG_F == KZ & P_CONF %in% PCONFS)
  CVs$genS_qLL <- filter(iStabi::data_CVs$genS_qLL, 
                         MOMENT_CONDITIONS == KZ & STRONGLY_IDENTIFIED_PARAMETERS == NSIPARAMS & P_CONF %in% PCONFS)
  CVs$S <- tibble('P_CONF' = PCONFS,
                 'CUTOFF_TOP' = map_dbl(PCONFS, ~qchisq(.x, KZMKX)))
  
  CVs
}
#CV$S <- tibble('SIG' = c(0.1, 0.05, 0.01, 0),
#               'VALUE' = c(0, qchisq(.9, KZMKX), 
#                           qchisq(.95, KZMKX), 
#                           qchisq(.99, KZMKX)))

CVs <- get_CVs()

Critical values are:

S

CVs$S

# A tibble: 2 × 2
  P_CONF CUTOFF_TOP
   <dbl>      <dbl>
1   0.85       2.07
2   0.9        2.71

qLLStab

CVs$qLLStab

# A tibble: 2 × 3
   DG_F P_CONF CUTOFF_TOP
  <dbl>  <dbl>      <dbl>
1     1   0.85       6.46
2     1   0.9        7.17

genS_qLL

CVs$genS_qLL

# A tibble: 2 × 4
  MOMENT_CONDITIONS STRONGLY_IDENTIFIED_PARAMETERS P_CONF CUTOFF_TOP
              <dbl>                          <dbl>  <dbl>      <dbl>
1                 1                              0   0.85       7.74
2                 1                              0   0.9        8.59

Data for plotting

build_set <- function(setName) {
  # setName = 'qLLStab'
  thisCV <- CVs[[setName]]
  statMax = filter(thisCV, P_CONF == PCONFS[2])$CUTOFF_TOP
  thisPconf = thisCV$P_CONF
  thisCV = thisCV$CUTOFF_TOP
  # 546 rows
  
  theseCols <- c(PARAMS_CFG$name, setName)
  aSet = Grid[get(setName) <= statMax, ..theseCols][, 
                                                    sigAt := thisPconf[ 1+findInterval(get(setName), thisCV,left.open=TRUE)]]
  aSet
  
}

Sets <- list('qLLStab' = build_set('qLLStab'),
             'S' = build_set('S'),
             'genS_qLL' = build_set('genS_qLL'))

PLOTTING

Grid for plotting

Force the starting values to be a multiple of the step size. This avoids problems and is clear what the true parameter space is.

Plots

True parameter vector $ =(= -0.75, = 0.15)$

Global search

Scatterplot for all searched values (useful for the global search)

S plot

# INTERACTIVE PLOTs ----
#### S-plot ----
pSbase <- plot_ly(data = Sets$S[sigAt<=0.90,], 
                 x=~alpha, 
                 y=~beta,
                 color='grey') %>%
  layout(title = 'S-set (90% confidence)',
         xaxis = list(title=plotly::TeX('\\alpha'), 
                        tickvals=seq(from=-1, to=1, by=0.20),
                        range = c(-1, 1),
                        fixedrange = TRUE),
           yaxis = list(title=plotly::TeX('\\beta'),  
                        tickvals=seq(from=-1, to=1, by=0.20),
                        range = c(-1, 1),
                        fixedrange = TRUE)) %>%
  add_trace(x = -0.75,y = 0.15,
    type = "scatter",   mode = "markers",
    marker = list(symbol = "square",  size = 12,  color = "red" ),
    inherit = FALSE,
    showlegend = FALSE) %>%
  config(mathjax = "cdn")

pSbase %>%
  add_trace(type='scatter', mode='markers', marker = list(color='blue'))

Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

qLL-Stab plot

# INTERACTIVE PLOTs ----
#### S-plot ----
pSbase <- plot_ly(data = Sets$qLLStab[sigAt<=0.90,], 
                 x=~alpha, 
                 y=~beta,
                 color='grey')%>%
  layout(title = 'qLLStab-set (90% confidence)',
         xaxis = list(title=plotly::TeX('\\alpha'), 
                        tickvals=seq(from=-1, to=1, by=0.20),
                        range = c(-1, 1),
                        fixedrange = TRUE),
           yaxis = list(title=plotly::TeX('\\beta'),  
                        tickvals=seq(from=-1, to=1, by=0.20),
                        range = c(-1, 1),
                        fixedrange = TRUE)) %>%
  add_trace(x = -0.75,y = 0.15,
    type = "scatter",   mode = "markers",
    marker = list(symbol = "square",  size = 12,  color = "red" ),
    inherit = FALSE,
    showlegend = FALSE) %>%
  config(mathjax = "cdn")

pSbase %>%
  add_trace(type='scatter', mode='markers', marker = list(color='blue'))

Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

pBase <- ggplot(data=Sets$qLLStab[sigAt<=0.90,], 
                mapping=aes(x=alpha, y=beta)) +
  scale_x_continuous(name=TeX(r"(\\alpha)"), limits=c(-1, 1), 
                     breaks=seq(from=-1, to=1, by=0.1)) +
  scale_y_continuous(name=TeX(r'($\\beta$)'), limits=c(-1, 1), 
                     breaks=seq(from=-1, to=1, by=0.1)) + 
  theme_bw() + theme(panel.grid=element_line(colour='#999999', linetype='14'), 
                     panel.grid.minor=element_line(colour='white'), 
                     text=element_text(size=10), 
                     axis.title.x = element_text(colour='red'),
                     axis.title.y = element_text(colour='blue', size=rel(0.9), hjust=0.9),
                     axis.text.x = element_text(size=rel(1.25), angle=90, vjust=0.5), 
                     axis.text.y = element_text(size=rel(1.25))) +
  geom_vline(xintercept=0, colour = 'black', linetype = 'solid') +
  geom_hline(yintercept=0, colour = 'black', linetype = 'solid') +   
  # Highlightinh the TRUE value
  geom_point(aes(x=-0.75, y=0.15), size=3, colour='black', fill='yellow')

# Scatterplot
pGplot <- pBase + geom_point(aes(colour=qLLStab))
ggsave('Unidentified_model_qLLStab_set.png', plot=pGplot, width=6, height=5)

Warning in geom_point(aes(x = -0.75, y = 0.15), size = 3, colour = "black", : All aesthetics have length 1, but the data has 546 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.

pGplot

Warning in geom_point(aes(x = -0.75, y = 0.15), size = 3, colour = "black", : All aesthetics have length 1, but the data has 546 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.